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Abstract 

Our paper computationally explores Ricker’s predator satiation model with the objective 
of studying how the extinction dynamics of an animal species having a two-stage life- 
cycle is affected by a sudden spike in mortality due to an extraneous extreme event. Our 
simulation model has been designed and implemented using sockeye salmon population 
data based on a stochastic version of Ricker’s model; with the shock size being reflected 
by a sudden reduction in the carrying capacity of the environment for this species. Our 
results show that even for a relatively marginal increase in the negative impact of an 
extreme event on the carrying capacity of the environment, a species with an otherwise 
stable population may be driven close to extinction. 
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Background and research objective 

PVA approaches do not normally consider the risk of catastrophic extreme events under 
the pretext that no population size can be large enough to guarantee survival of a species 
in the event of a large-scale natural catastrophe. '■'^Nevertheless, it is only very intuitive 
that some species are more “delicate” than others; and although not presently under any 
clearly observed threat, could become threatened with extinction very quickly if an 
extreme event was to occur even on a low-to-moderate scale. The term “extreme event” is 
preferred to “catastrophe” because catastrophe usually implies a natural event whereas; 
quite clearly; the chance of man-caused extreme events poses a much greater threat at 
present to a number of animal species as compared to any large-scale natural catastrophe. 



An animal has a two-stage life eyele when; in the first stage, newborns beeome immature 
youths and in the seeond stage; the immature youths beeome mature adults. Therefore, in 
terms of the stage-speeifie approaeh, if Yt denotes the number of immature young in 
stage t and At denotes the number of mature adults, then the number of adults in year t + 
1 will be some proportion of the young, speoifieally those that survive to the next 
(reproduetive) stage. Then the formal relationship between the number of mature adults 
in the next stage and the number of immature youths at present may be written as 
follows: 



At+i=XYt 

Here A, is the survival probability, i.e. it is the probability of survival of a youth to 
maturity. The number of young next year will depend on the number of adults in t: 

Yt+i = f(At) 

Here f deseribes the reproduetion relation between mature adults and next year’s young. 

This is a straightforward system of simultaneous differenee equations whieh may be 
analytieally solved using a variation of the cobwebbing approach. The solution proeess 
begins with an initial point (Yi, Ai) and iteratively determines the next point (Y2, A2). If 
predator satiation is built into the proeess, then we simply end up with Rieker’s model: 

'\r A — At/K 

A t + 1 ~ ot AfC 

Here a is the maximum reproduetion rate (for an initial small population) and K is the 
population size at whieh the reproduetion rate is approximately half its maximum 
Putting fi = Vk we ean re-write Rieker’s equation as follows: 

Yt+i =aAte"^^‘ 



It has been shown that if (Yo, Ao) lies within the first of three possible ranges, (Yn, A„) 
approaehes (0, 0) in sueeessive years and the population beeomes extinet. If (Yo, Ao) lies 
within the third range then (Yn, A„) equilibrate to a steady-state value of (Y*, A*). 
Populations that begin with (Yo, Aq) within the seeond range oseillate between (Y*, 0) 
and (0, A*). Sueh alternating behavior indieates one of the year elasses, or eohorts, 
beeome extinet while the other persists i.e. adult breeding stoek appear only every other 
year. Thus the model reveals that three quite different results oeeur depending initially 
only on the starting sizes of the population and its distribution among the two stages. 

We use the same basie model in our researeh but instead of analytieally solving the 
system of differenee equations, we use the same to simulate the population dynamies as a 
stoehastie proeess implemented on an MS-Exeel spreadsheet. Rather than using a elosed- 
form equation like Rieker’s model to represent the funetional relationship between Yt+ 1 
and At, we use a Monte Carlo method to simulate the stage-transition proeess within 




Ricker’s framework; introdueing a massive perturbation with a very small probability in 
order to emulate a eatastrophie event 



Conceptual framework 

We have a formulated a stoehastic population growth model with an inbuilt eapaeity to 
generate an extreme event based on a theoretieal probability distribution. The non- 
stochastic part of the model eorresponds to Ricker’s relationship between Yt + i and At. 
The stoehastie part has to do with whether or not an extreme event oeeurs at a particular 
time point. The gamma distribution has been chosen to make the probability distribution 
for the extreme event a skewed one as it is likely to be in reality. Instead of analytically 
solving the system of simultaneous differenee equations iteratively in some variation of 
the eobwebbing method, we have used them in a spreadsheet model to simulate the 
population growth over a span of ten time periods. 

We apply a eomputational methodology whereby the initial number of immature young is 
hypothesized to either attain the expeeted number predieted by Rieker’s model or 
drastieally fall below that number at the end of every stage, depending on whether an 
extraneous extreme event does not oecur or actually occurs. The mortalities as a result of 
an extreme event at any time point is expressed as a pereentage of the pristine population 
size for a elearer eomparative view. 



Model building 

Among various faunal speeies, the population dynamies of the soekeye salmon 
(oncorhynchus nerka) has been most extensively studied using Riekert’s model. Salmon 
are unique in that they breed in partieular fresh water systems before they die. Their 
offspring migrates to the oeean and upon reproduetive maturity, they are guided by a 
hitherto unaeeounted instinctive drive to swim baek to the very same fresh waters where 
they were bom to spawn their own offspring and perish. Salmon populations thus are 
very sensitive to habitat ehanges and human aetivities that have a negative impaet on 
riparian eeosystems that serve as breeding grounds for salmon ean adversely affeet the 
peculiar life-cycle of the salmon. Many of the aneient salmon mns (notably those in 
California river systems) have now gone extinct and it is our hypothesis that an even 
seemingly stabilized population ean be rapidly driven to extinetion due to the effeet of an 
extraneous (quite possibly man-made) extreme event with the eapaeity to eause mass 
mortality. The following table shows the four-year averages of the soekeye salmon 
population in the Skeena river system in British Columbia in the first half of the twentieth 
eentury. 




Year 


Population (in thousands) 


1908 


1,098 


1912 


740 


1916 


714 


1920 


615 


1924 


706 


1928 


510 


1932 


278 


1936 


448 


1940 


528 


1944 


639 


1948 


523 



{Source: http://www- 

rohan.sdsu.edu/~jmahaJJy/courses/sOO/mathl 2 l/lectures/product_rule/product.html#Rick 

er'sModel) 



A non-linear least squares best-fit to Ricker’s model is obtained for the above set of data 
is obtained as follows: 



n 

Minimize , where dt is the actual population size in year t. 

(=1 



The necessary conditions to the above least squares best-fit problem is obtained as 
follows: 



5(s^)/5a = Oi whereby we get a* ~ 1.54 and fl* ~ 7.8 x 10“^ 

Plugging these parameters into Ricker’s model indeed yields a fairly good approximation 
of the salmon population stabilization in the Skeena river system in the first half of the 
previous century. 

As the probability distribution of an extraneous extreme event is likely to be a highly 
skewed one, we have generated our random variables from the cumulative distribution 
function {cdj) of the gamma distribution rather than the normal distribution. The 
distribution boundaries are fixed by generating random integers in the range 1 to 1 00 and 
using these random integers to define the shape and scale parameters of the gamma 
distribution. The gamma distribution performs better than the normal distribution when 
the distribution to be matched is highly right-skewed; as is desired in our model. The 
combination of a large variance and a lower limit at zero makes fitting a normal 
distribution rather unsuitable in such cases. The probability density function of the 
gamma distribution is given as follows: 





f (x, a, b) = {b“T{a)} 'x“ 'e for x > 0 



Here a > 0 is the shape parameter and P > 0 is the seale parameter of the gamma 
distribution. The eumulative distribution function may be expressed in terms of the 
incomplete gamma function as follows: 



F (x, a, b) = j / {u)du = /(a, x/b)/ F(a) 

0 

In our spreadsheet model, we have F {R, R/2, 2) as our cdf of the gamma distribution. 
Here R is an integer randomly sampled from the range 1 to 100. An interesting statistical 
result of having these values for x, a and P is that the cumulative gamma distribution 
value becomes equalized with the value [1 - x (^)] having R degrees of freedom, thus 
allowing goodness-offit tests. 

Our model is specifically designed to simulate the extinction dynamics of sockeye 
salmon population using a stochastic version of Ricker’s model; with the shock size 
being based on a sudden reduction in the parameter K i.e. the carrying capacity of the 
environment for this species. The model parameters are same as those of Ricker’s model 
i.e. a and P (which is the reciprocal of K). We have kept a constant at all times at 1.54, 
which was the least squares best-fit value obtained for that parameter. We have kept a P 
of 0.00078 (i.e. the best-lit value) when no extreme event occurs and have varied the P 
between 0.00195 and 0.0156 (i.e. between 2.5 times to 20 times the best- fit value) for 
cases where an extreme event occurred. We have a third parameter c which is basically a 
‘switching constant’ that determines whether an extreme event occurs or not. The switch 
is turned on triggering an extreme event when a random draw from a cumulative gamma 
distribution yields a value less than or equal to c. Using F {R, R/2, 2) as our cdf of the 
gamma distribution where 7? is a randomly drawn integer in the range (1, 100) means that 
the cumulative gamma function will randomly select from the approximate interval 0.518 
~ 0.683. By fixing the value of c at 0.5189 in our model we have effectively reduced the 
probability of occurrence of an extreme event to a miniscule magnitude relative to that of 
an extreme event not occurring. We have used the sockeye salmon population data from 
the table presented earlier For each level of the P parameter, we simulated the system and 
observed the maximum possible number of mortalities from an extreme event at that 
level of p. The results are reported below. 



Results obtained from the simulation model 




Worst-case effect of extreme event on sockeye salmon population 
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We made 100 independent simulation runs for each of the eight levels of p. The low 
probability of extreme event assigned in our study yielded a mean of 1.375 for the 
number of observed worst-case scenarios (i.e. situations of maximum mortality) with a 
standard deviation of approximately 0.92. The worst-case scenarios for our choice of 
parameters necessarily occur if the extreme event occurs in the first time point when the 
species population is at its maximum size. Our model shows that in worst-case scenarios, 
the size of surviving population after an extreme event that could seed the ultimate 
recovery of the species to pre-catastrophe numbers (staying within the broad framework 
of Ricker’s model) drops from about 18% of the pristine population size for a shock size 
corresponding to 2.5 times the best-fit P; to only about 0.000005% of the pristine 
population size for a shock size corresponding to 20 times the best-fit p. 

Therefore, if the minimum required size of the surviving population is at least say 20% of 
the pristine population in order to survive and recover to pre-catastrophe numbers, the 
species could go extinct if an extreme event caused a little more than two-fold decrease in 
the environmental carrying capacity! Even if the minimum required size for recovery was 
relatively low at say around 2% of the pristine population, an extreme event that caused a 
five-fold decrease in the environmental carrying capacity could very easily force the 
species to the brink of extinction. An immediate course of future extension of our work 
would be allowing the fecundity parameter a to be affected by extreme events as is very 
likely in case of say a large-scale chemical contamination of an ecosystem due to a faulty 
industrial waste-treatment facility. 



Conclusion 



Our study has shown that even for a relatively marginal 2.5-fold deerease in the 
environmental carrying capacity due to an extreme event, a worst-case scenario could 
mean a mortality figure well above 80% of the pristine population. As a guide for future 
PVA studies we may suggest that one should not be deterred simply by the notion that 
extreme events are uncontrollable and hence outside the purview of computational 
modeling. Indeed the effect of an extreme event can almost always prove to be fatal for a 
species but nevertheless, as our study shows, there is ample scope and justification for 
future scientific enquiries into the relationship between survival probability of a species 
and the adverse impact of an extreme event on ecological sustainability. 
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